---
title: "Van Egeren et al (2018) Supplementary Information"
output: html_notebook
---

```{r echo=FALSE}
library(readr)
library(ggplot2)
library(reshape2)
library(RColorBrewer)
source("./utils.R")
```

## Supplementary Fig. 1a

```{r include=FALSE}
filepath = "./"
all_dat = c()
lifetimes = c()
input_vars = c()
for (k in 2:6){
  for (i in c("heritable")){
    lifetimes = c(lifetimes, i)
    input_vars = c(input_vars, k)
    for (j in 1:1000){
      filename = paste(filepath, "lognorm_ss_dist/", i, "/", k, "/fit_sim_",j,"type_0.oevo", sep="")
      temp = read_csv(filename, skip=1, col_names=FALSE)
      est_mean = mean(unlist(t(temp[1,2:101])))
      to_deliver = unlist(t(temp[1,2:101]))/est_mean
      all_dat = c(all_dat, to_deliver)
    }
  }
}
dist_data = data.frame(fit=all_dat, lifetimes=rep(lifetimes,each=1000), input_variances=rep(input_vars,each=100000))
```

```{r echo=FALSE}
g1 = plot_dists_var(dist_data)
g1 + xlim(0.5, 1.5)
```

## Supplementary Fig. 1b

```{r include=FALSE}
filepath = "./"
# THIS WILL LIKELY TAKE A WHILE
variances = extract_vars(filepath)
her_variances = extract_vars_heritable(filepath)
```

```{r echo=FALSE}
ggplot(data=variances) +geom_point(aes(x=-input_var,y=log_var,color=lifetimes)) + geom_errorbar(aes(x=-input_var, ymin=log_lower, ymax=log_upper,color=lifetimes), width=0.05) + theme_bw() + theme(legend.position = "none", axis.title.x=element_blank(), axis.title.y=element_blank(), text = element_text(size=20), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"), panel.border = element_blank()) + scale_color_gradient2(low="#cccc00", mid="#ff3399", high="#9900ff", midpoint=4)  + geom_point(data=her_variances, aes(x=-her_variances$input_var,y=her_variances$log_var), color="#4400cc") + geom_errorbar(data=her_variances, aes(x=-her_variances$input_var, ymin=log_lower, ymax=log_upper), color="#4400cc", width=0.05)
```

## Supplementary Fig. 2a-c

```{r include=FALSE}
filepath = "./"
all_dat_twopt = c()
lifetimes = c()
input_vars = c()
for (k in 4:4){
  for (i in c(1,3,5,10,50,"heritable")){
    lifetimes = c(lifetimes, i)
    input_vars = c(input_vars, k)
    for (j in 1:1000){
      filename = paste(filepath, "data_twopt_regress_dist/", i, "/", k, "/fit_sim_",j,"type_0.oevo", sep="")
      temp = read_csv(filename, skip=1, col_names=FALSE)
      est_mean = mean(unlist(t(temp[1,2:101])))
      to_deliver = unlist(t(temp[1,2:101]))/est_mean
      all_dat_twopt = c(all_dat_twopt, to_deliver)
    }
  }
}
twopt_data = data.frame(fit=all_dat_twopt, lifetimes=rep(lifetimes,each=100000), input_variances=rep(input_vars,each=1000))
```

```{r include=FALSE}
all_dat_gamexp = c()
lifetimes = c()
shape = c()
for (k in c("gamma", "expo")){
  for (i in c(1,3,5,10,50,"heritable")){
    lifetimes = c(lifetimes, i)
    shape = c(shape, k)
    for (j in 1:1000){
      filename = paste(filepath, "data_gam_exp_regress_dist/", k, "/", i, "/fit_sim_",j,"type_0.oevo", sep="")
      temp = read_csv(filename, skip=1, col_names=FALSE)
      est_mean = mean(unlist(t(temp[1,2:101])))
      to_deliver = unlist(t(temp[1,2:101]))/est_mean
      all_dat_gamexp = c(all_dat_gamexp, to_deliver)
    }
  }
}
gamexp_data = data.frame(fit=all_dat_gamexp, lifetimes=rep(lifetimes,each=1000), input_dist=rep(shape,each=100000))
gamexp_data$input_dist = c(rep("gamma", times=length(gamexp_data$input_dist)/2),rep("expo", times=length(gamexp_data$input_dist)/2))
gamexp_data$lifetimes = rep(rep(c(1,3,5,10,50,"heritable"), each=100000), 2)
```

```{r echo=FALSE}
# Supplementary Fig. 2a
dist_data = gamexp_data[gamexp_data$input_dist == "gamma",]
g1 = plot_dists(dist_data)
g1

# Supplementary Fig. 2b
dist_data = gamexp_data[gamexp_data$input_dist == "expo",]
g1 = plot_dists(dist_data)
g1

# Supplementary Fig. 2c
dist_data = twopt_data
g1 = plot_dists(dist_data)
g1
```

## Supplementary Fig. 2d

```{r include=FALSE}
filepath = "./"
all_dat_lgnorm = c()
lifetimes = c()
input_vars = c()
for (k in 4:4){
  for (i in c("heritable")){
    lifetimes = c(lifetimes, i)
    input_vars = c(input_vars, k)
    for (j in 1:1000){
      filename = paste(filepath, "lognorm_ss_dist/", i, "/", k, "/fit_sim_",j,"type_0.oevo", sep="")
      temp = read_csv(filename, skip=1, col_names=FALSE)
      est_mean = mean(unlist(t(temp[1,2:101])))
      to_deliver = unlist(t(temp[1,2:101]))/est_mean
      all_dat_lgnorm = c(all_dat_lgnorm, to_deliver)
    }
  }
}
lgnorm_data = data.frame(fit=all_dat_lgnorm, lifetimes=rep(lifetimes,each=100000), input_variances=rep(input_vars,each=1000))
```

```{r echo=FALSE}
ggplot() + geom_line(aes(x=gamexp_data[gamexp_data$input_dist == "gamma" & gamexp_data$lifetimes == "heritable",]$fit), stat="density", color="orange", size=0.8) + geom_line(aes(x=gamexp_data[gamexp_data$input_dist == "expo" & gamexp_data$lifetimes == "heritable",]$fit), stat="density", color="purple", size=0.8) + theme_bw() + geom_line(aes(x=twopt_data[twopt_data$lifetimes == "heritable",]$fit), stat="density", color="green3", size=0.8) + theme(legend.position = "none", axis.title.x=element_blank(), axis.title.y=element_blank(), text = element_text(size=20), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"), panel.border = element_blank())+ geom_line(aes(x=lgnorm_data[lgnorm_data$lifetimes == "heritable",]$fit), stat="density", color="black", size=0.8)
```

## Supplementary Fig. 3a-b

```{r include=FALSE}
filepath = "./data_r1_regress_twopt/"
first = c(0.9, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1)
second = c(1,5,10,50,100)
dat = extract_data(filepath, first, second, 10000)
```

```{r echo=FALSE}
plots = plot_fix_tunnel(dat$fix_dat, dat$simple_dat, dat$her_dat)
# Supplementary Fig. 3a
plots$fix

# Supplementary Fig. 3b
plots$tunnel
```

## Supplementary Fig. 3c-d

```{r include=FALSE}
filepath = "./data_r1_regress_exp/"
first = c(0.9, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1)
second = c(1,5,10,50,100)
dat = extract_data(filepath, first, second, 10000)
```

```{r echo=FALSE}
plots = plot_fix_tunnel(dat$fix_dat, dat$simple_dat, dat$her_dat)
# Supplementary Fig. 3c
plots$fix

# Supplementary Fig. 3d
plots$tunnel
```

## Supplementary Fig. 4a-b

```{r include=FALSE}
filepath = "./data_u1_regress/"
first = c(2,3,4,5,6)
second = c(1,5,10,50,100)
dat = extract_data(filepath, first, second, 10000)
dat$fix_dat$first = -dat$fix_dat$first
dat$simple_dat$first = -dat$simple_dat$first
dat$her_dat$first = -dat$her_dat$first
```

```{r echo=FALSE}
plots = plot_fix_tunnel(dat$fix_dat, dat$simple_dat, dat$her_dat)
# Supplementary Fig. 4a
plots$fix + scale_x_continuous(breaks=c(-6, -5, -4, -3, -2), labels=c("1e-6", "1e-5", "1e-4", "1e-3", "1e-2"))

# Supplementary Fig. 4b
plots$tunnel + scale_x_continuous(breaks=c(-6, -5, -4, -3, -2), labels=c("1e-6", "1e-5", "1e-4", "1e-3", "1e-2"))
```

## Supplementary Fig. 4c-d

```{r include=FALSE}
filepath = "./data_u2_regress/"
first = c(2,3,4,5,6)
second = c(1,5,10,50,100)
dat = extract_data(filepath, first, second, 10000)
dat$fix_dat$first = -dat$fix_dat$first
dat$simple_dat$first = -dat$simple_dat$first
dat$her_dat$first = -dat$her_dat$first
```

```{r echo=FALSE}
plots = plot_fix_tunnel(dat$fix_dat, dat$simple_dat, dat$her_dat)
# Supplementary Fig. 4c
plots$fix + scale_x_continuous(breaks=c(-6, -5, -4, -3, -2), labels=c("2e-6", "2e-5", "2e-4", "2e-3", "2e-2"))

# Supplementary Fig. 4d
plots$tunnel + scale_x_continuous(breaks=c(-6, -5, -4, -3, -2), labels=c("2e-6", "2e-5", "2e-4", "2e-3", "2e-2"))
```

## Supplementary Fig. 5a-b

```{r include=FALSE}
filepath = "./data_r2_regress/"
first = c(1, 1.1, 1.5, 2, 3, 4, 5)
second = c(1,5,10,50,100)
dat = extract_data(filepath, first, second, 10000)
```

```{r echo=FALSE}
plots = plot_fix_tunnel(dat$fix_dat, dat$simple_dat, dat$her_dat)
# Supplementary Fig. 5a
plots$fix

# Supplementary Fig. 5b
plots$tunnel
```

## Supplementary Fig. 5c

```{r include=FALSE}
filepath = "./"
dat3 = data.frame(A=numeric(0), B=numeric(0), C=numeric(0), D=numeric(0))
for (i in c(0.8, 0.9, 0.95, 0.99, 1)){
  for (j in c(1, 1.01, 1.05, 1.1, 1.5, 2, 5)){
    filename = paste(filepath, "data_r2_r1/", i, "/", j, "/heritable/type_1_tunnel.oevo", sep="")
    temp = read_csv(filename, col_names=FALSE)
    tunnel_prob = mean(temp$X2)
    filename = paste(filepath, "data_r2_r1/", i, "/", j, "/heritable/iftype2.oevo", sep="")
    temp = read_csv(filename, col_names=FALSE)
    fix_prob = mean(temp$X2)
    dat3 = rbind(dat3, data.frame(A=i, B=tunnel_prob, C=fix_prob, D=j))
  }
}
colnames(dat3) <- c("r1", "tunnel", "fix", "r2")
dat3["tunnel_norm"] = dat3$tunnel/dat3$fix
```

```{r include=FALSE}
filepath = "./"
dat4 = data.frame(A=numeric(0), B=numeric(0), C=numeric(0), D=numeric(0))
for (i in c(0.8, 0.9, 0.95, 0.99, 1)){
  for (j in c(1, 1.01, 1.05, 1.1, 1.5, 2, 5)){
    filename = paste(filepath, "data_r2_r1/", i, "/", j, "/simple/type_1_tunnel.oevo", sep="")
    temp = read_csv(filename, col_names=FALSE)
    tunnel_prob = mean(temp$X2)
    filename = paste(filepath, "data_r2_r1/", i, "/", j, "/simple/iftype2.oevo", sep="")
    temp = read_csv(filename, col_names=FALSE)
    fix_prob = mean(temp$X2)
    dat4 = rbind(dat4, data.frame(A=i, B=tunnel_prob, C=fix_prob, D=j))
  }
}
colnames(dat4) <- c("r1", "tunnel", "fix", "r2")
dat4["tunnel_norm"] = dat4$tunnel/dat4$fix
```

```{r include=FALSE}
filepath = "./"
dat5 = data.frame(A=numeric(0), B=numeric(0), C=numeric(0), D=numeric(0))
for (i in c(0.8, 0.9, 0.95, 0.99, 1)){
  for (j in c(1, 1.01, 1.05, 1.1, 1.5, 2, 5)){
    filename = paste(filepath, "data_r2_r1/", i, "/", j, "/type/type_1_tunnel.oevo", sep="")
    temp = read_csv(filename, col_names=FALSE)
    tunnel_prob = mean(temp$X2)
    filename = paste(filepath, "data_r2_r1/", i, "/", j, "/type/iftype2.oevo", sep="")
    temp = read_csv(filename, col_names=FALSE)
    fix_prob = mean(temp$X2)
    dat5 = rbind(dat5, data.frame(A=i, B=tunnel_prob, C=fix_prob, D=j))
  }
}
colnames(dat5) <- c("r1", "tunnel", "fix", "r2")
dat5["tunnel_norm"] = dat5$tunnel/dat5$fix
```

```{r echo=FALSE}
dat3["difference"] = dat3$tunnel_norm - dat4$tunnel_norm
to_plot = dat3[dat3$r2 < 4,]

smoothed = loess (difference ~ r2 * r1, to_plot, degree = 0, span = 0.1)
x = (seq(min((to_plot$r2)), max((to_plot$r2)), 0.005))
y = seq(min(to_plot$r1), max(to_plot$r1), 0.005)
interpolated <- predict (smoothed, expand.grid (r2 =x, r1 = y))
melted = melt(interpolated)
melted$r2 = x
melted$r1 = rep(y, times=rep(length(x), length(y)))
ggplot() + geom_tile(data=melted, aes(x=(r2), y=r1, fill=value)) + scale_fill_gradient2(limits=c(-0.5, 0.1), low="red", high="green") + theme_minimal()  + ylim(min(to_plot$r1), max(to_plot$r1)) + geom_point(data=to_plot, aes(x=(r2), y=r1), size= 1.7) + geom_point(data=to_plot, aes(x=(r2), y=r1, color=difference), size=1.3) + scale_color_gradient2(limits=c(-0.5, 0.1), low="red", high="green")+ theme(panel.grid = element_blank(), panel.border = element_blank())  + theme(text = element_text(size=20))+ xlab("final stage fitness r2") + ylab("valley depth r1")
```

```{r echo=FALSE}
dat5["difference"] = dat5$tunnel_norm - dat4$tunnel_norm
to_plot = dat5[dat5$r2 < 4,]

smoothed = loess (difference ~ r2 * r1, to_plot, degree = 0, span = 0.1)
x = (seq(min((to_plot$r2)), max((to_plot$r2)), 0.005))
y = seq(min(to_plot$r1), max(to_plot$r1), 0.005)
interpolated <- predict (smoothed, expand.grid (r2 =x, r1 = y))
melted = melt(interpolated)
melted$r2 = x
melted$r1 = rep(y, times=rep(length(x), length(y)))
ggplot() + geom_tile(data=melted, aes(x=(r2), y=r1, fill=value)) + scale_fill_gradient2(limits=c(-0.5, 0.1), low="red", high="green") + theme_minimal()  + ylim(min(to_plot$r1), max(to_plot$r1)) + geom_point(data=to_plot, aes(x=(r2), y=r1), size= 1.7) + geom_point(data=to_plot, aes(x=(r2), y=r1, color=difference), size=1.3) + scale_color_gradient2(limits=c(-0.5, 0.1), low="red", high="green")+ theme(panel.grid = element_blank(), panel.border = element_blank())  + theme(text = element_text(size=20))+ xlab("final stage fitness r2") + ylab("valley depth r1")
```

## Supplementary Fig. 6

```{r include=FALSE}
filepath = "./data_N_regress/"
first = c(10, 20, 50, 100, 200, 500)
second = c(1,5,10,50,100)
dat = extract_data(filepath, first, second, 10000)
```

```{r echo=FALSE}
plots = plot_fix_tunnel(dat$fix_dat, dat$simple_dat, dat$her_dat)
# Supplementary Fig. 6a
plots$fix

# Supplementary Fig. 6b
plots$tunnel
```

## Supplementary Fig. 7

```{r include=FALSE}
filepath = "./data_r1_regress_stoch/"
first = c(0.9, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1)
second = c(1,3,5,10,50)
dat = extract_data(filepath, first, second, 10000)
```

```{r echo=FALSE}
fix_dat = dat$fix_dat
her_dat = dat$her_dat
simple_dat = dat$simple_dat

# Supplementary Fig. 7a
ggplot(data=fix_dat) + geom_line(aes(x=first, y=fix, color=sapply(second, FUN=to_colorscale2), group=second)) + theme_bw() + geom_errorbar(aes(x=first,ymin=fix-fix_SEM*1.96, ymax=fix+fix_SEM*1.96, color=sapply(second, FUN=to_colorscale2)), width=0) + geom_line(data=simple_dat, aes(x=first, y=fix), color="#707070", size=1) + geom_errorbar(data=simple_dat, aes(x=first, ymin=fix-fix_SEM*1.96, ymax=fix+fix_SEM*1.96), color="#707070", width=0) + geom_line(data=her_dat, aes(x=first, y=fix), color="#4400cc", size=1) + geom_errorbar(data=her_dat, aes(x=first, ymin=fix-fix_SEM*1.96, ymax=fix+fix_SEM*1.96), color="#4400cc", width=0) + scale_color_gradient2(low="#cccc00", mid="#ff3399", high="#9900ff", midpoint=4) + theme(legend.position = "none", axis.title.x=element_blank(), axis.title.y=element_blank(), text = element_text(size=20), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"), panel.border = element_blank())

# Supplementary Fig. 7b
ggplot(data=fix_dat) + geom_line(aes(x=first, y=tunnel, color=sapply(second, FUN=to_colorscale2), group=second)) + theme_bw() + geom_errorbar(aes(x=first,ymin=tunnel-tunnel_SEM*1.96, ymax=tunnel+tunnel_SEM*1.96, color=sapply(second, FUN=to_colorscale2)), width=0) + geom_line(data=simple_dat, aes(x=first, y=tunnel), color="#707070", size=1) + geom_errorbar(data=simple_dat, aes(x=first, ymin=tunnel-tunnel_SEM*1.96, ymax=tunnel+tunnel_SEM*1.96), color="#707070", width=0) + geom_line(data=her_dat, aes(x=first, y=tunnel), color="#4400cc", size=1) + geom_errorbar(data=her_dat, aes(x=first, ymin=tunnel-tunnel_SEM*1.96, ymax=tunnel+tunnel_SEM*1.96), color="#4400cc", width=0)+ scale_color_gradient2(low="#cccc00", mid="#ff3399", high="#9900ff", midpoint=4) + theme(legend.position = "none", axis.title.x=element_blank(), axis.title.y=element_blank(), text = element_text(size=20), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"), panel.border = element_blank())
```

## Supplementary Fig. 8

```{r include=FALSE}
filepath = "./data_rev/"
her_variances = extract_vars_heritable(filepath)

dat = data.frame(A=numeric(0), B=numeric(0), C=numeric(0), D=numeric(0))
for (i in 2:6){
  for (j in c(0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99)){
    filename = paste(filepath, "data_var_r1/", j, "/", i, "/heritable/type_1_tunnel.oevo", sep="")
    temp = read_csv(filename, col_names=FALSE)
    tunnel_prob = mean(temp$X2)
    filename = paste(filepath, "data_var_r1/", j, "/", i, "/heritable/iftype2.oevo", sep="")
    temp = read_csv(filename, col_names=FALSE)
    fix_prob = mean(temp$X2)
    dat = rbind(dat, data.frame(A=i, B=tunnel_prob, C=fix_prob, D=j))
  }
}
colnames(dat) <- c("var", "tunnel", "fix", "r1")
dat["tunnel_norm"] = dat$tunnel/dat$fix
dat["var_est"] = rep(her_variances$mean_var, times=c(9,9,9,9,9))
dat["log_lower"] = rep(her_variances$log_lower, times=c(9,9,9,9,9))
dat["log_var"] = rep(her_variances$log_var, times=c(9,9,9,9,9))
dat["log_upper"] = rep(her_variances$log_upper, times=c(9,9,9,9,9))

```
```{r echo=FALSE}
# Supplementary Fig. 8a
smoothed = loess (tunnel_norm ~ log_var * r1, dat, degree = 0, span = 0.1)
x = (seq(min(dat$log_var), max(dat$log_var), 0.005))
y = seq(min(dat$r1), max(dat$r1), 0.005)
interpolated <- predict (smoothed, expand.grid (log_var =x, r1 = y))
melted = melt(interpolated)
melted$log_var = x
melted$r1 = rep(y, times=rep(length(x), length(y)))
ggplot() + geom_tile(data=melted, aes(x=log_var, y=r1, fill=value)) + scale_fill_gradientn(colours=rev(brewer.pal(9,"RdPu")), limits=c(0.18, 1.01)) + theme_minimal()  + ylim(min(dat$r1), max(dat$r1)) + geom_point(data=dat, aes(x=log_var, y=r1), size= 1.7) + geom_point(data=dat, aes(x=log_var, y=r1, color=tunnel_norm), size=1.3) + scale_color_gradientn(colours=rev(brewer.pal(9,"RdPu")), limits=c(0.18, 1.01))+ theme(panel.grid = element_blank(), panel.border = element_blank()) + scale_x_continuous(breaks = c(her_variances$log_var), labels=formatC(10^her_variances$log_var, format = "e", digits = 2), limits=c((min(dat$log_var)), (max(dat$log_var)))) + theme(text = element_text(size=20))+ xlab("avg fitness variance (log scale)") + ylab("valley depth r1")

```

```{r echo=FALSE}
# Supplementary Fig. 8b
smoothed = loess (fix ~ log_var * r1, dat, degree = 0, span = 0.1)
x = (seq(min(dat$log_var), max(dat$log_var), 0.005))
y = seq(min(dat$r1), max(dat$r1), 0.005)
interpolated <- predict (smoothed, expand.grid (log_var =x, r1 = y))
melted = melt(interpolated)
melted$log_var = x
melted$r1 = rep(y, times=rep(length(x), length(y)))
ggplot() + geom_tile(data=melted, aes(x=log_var, y=r1, fill=value)) + scale_fill_gradientn(colours=rev(brewer.pal(9,"RdPu")), limits=c(0, 0.5)) + theme_minimal()  + ylim(min(dat$r1), max(dat$r1)) + geom_point(data=dat, aes(x=log_var, y=r1), size= 1.7) + geom_point(data=dat, aes(x=log_var, y=r1, color=fix), size=1.3) + scale_color_gradientn(colours=rev(brewer.pal(9,"RdPu")), limits=c(0, 0.5))+ theme(panel.grid = element_blank(), panel.border = element_blank()) + scale_x_continuous(breaks = c(her_variances$log_var), labels=formatC(10^her_variances$log_var, format = "e", digits = 2), limits=c((min(dat$log_var)), (max(dat$log_var)))) + theme(text = element_text(size=20))+ xlab("avg fitness variance (log scale)") + ylab("valley depth r1")

```

```{r include=FALSE}
filepath = "./data_rev/"
dat2 = data.frame(A=numeric(0), B=numeric(0), C=numeric(0), D=numeric(0))
for (i in 1:6){
  for (j in c(0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99)){
    filename = paste(filepath, "data_var_r1/", j, "/", i, "/type/type_1_tunnel.oevo", sep="")
    temp = read_csv(filename, col_names=FALSE)
    tunnel_prob = mean(temp$X2)
    filename = paste(filepath, "data_var_r1/", j, "/", i, "/type/iftype2.oevo", sep="")
    temp = read_csv(filename, col_names=FALSE)
    fix_prob = mean(temp$X2)
    dat2 = rbind(dat2, data.frame(A=i, B=tunnel_prob, C=fix_prob, D=j))
  }
}
colnames(dat2) <- c("var", "tunnel", "fix", "r1")
dat2["tunnel_norm"] = dat2$tunnel/dat2$fix
dat2["sd_est"] = sqrt(10^-dat2$var)
dat2["var_est"] = 10^-dat2$var
```

```{r include=FALSE}
filepath = "./data_rev/"
dat_simple = data.frame(B=numeric(0), C=numeric(0), D=numeric(0))

  for (j in c(0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99)){
    filename = paste(filepath, "data_var_r1/", j, "/simple/type_1_tunnel.oevo", sep="")
    temp = read_csv(filename, col_names=FALSE)
    tunnel_prob = mean(temp$X2)
    filename = paste(filepath, "data_var_r1/", j, "/simple/iftype2.oevo", sep="")
    temp = read_csv(filename, col_names=FALSE)
    fix_prob = mean(temp$X2)
    dat_simple = rbind(dat_simple, data.frame(B=tunnel_prob, C=fix_prob, D=j))
  }
colnames(dat_simple) <- c("tunnel", "fix", "r1")
dat_simple["tunnel_norm"] = dat_simple$tunnel/dat_simple$fix
dat_simple["tunnel_SEM"] = sqrt((dat_simple$tunnel_norm) * (1-dat_simple$tunnel_norm) / (dat_simple$fix * 10000))
dat_simple["fix_SEM"] = sqrt((dat_simple$fix) * (1-dat_simple$fix) / (10000))
```

```{r echo=FALSE}
# Supplementary Fig. 8c
smoothed = loess (tunnel_norm ~ log10(var_est) * r1, dat2, degree = 0, span = 0.07)
x = 10^(seq(log10(min(dat2$var_est)), log10(max(dat2$var_est)), 0.005))
y = seq(min(dat2$r1), max(dat2$r1), 0.005)
x_tighter = seq(min(dat2$var_est), max(dat2$var_est), 0.01)
interpolated <- predict (smoothed, expand.grid (var_est =x, r1 = y))
melted = melt(interpolated)
melted$var_est = x
melted$r1 = rep(y, times=rep(length(x), length(y)))
ggplot() + geom_tile(data=melted, aes(x=log10(var_est), y=r1, fill=value)) + scale_fill_gradientn(colours=rev(brewer.pal(9,"RdPu")), limits=c(0.18, 1.01)) + theme_minimal()  + ylim(min(dat2$r1), max(dat2$r1)) + geom_point(data=dat2, aes(x=log10(var_est), y=r1), size= 1.7) + geom_point(data=dat2, aes(x=log10(var_est), y=r1, color=tunnel_norm), size=1.3) + scale_color_gradientn(colours=rev(brewer.pal(9,"RdPu")), limits=c(0.18, 1.01))+ theme(panel.grid = element_blank(), panel.border = element_blank()) + scale_x_continuous(breaks = c(unique(log10(dat2$var_est))), labels=formatC(c(unique((dat2$var_est))), format = "e", digits=0), limits=c(-6, -1))+ theme(text = element_text(size=20))+ xlab("fitness variance (log scale)") + ylab("valley depth r1")

```

```{r echo=FALSE}
# Supplementary Fig. 8d
smoothed = loess (fix ~ log10(var_est) * r1, dat2, degree = 0, span = 0.07)
x = 10^(seq(log10(min(dat2$var_est)), log10(max(dat2$var_est)), 0.005))
y = seq(min(dat2$r1), max(dat2$r1), 0.005)
x_tighter = seq(min(dat2$var_est), max(dat2$var_est), 0.01)
interpolated <- predict (smoothed, expand.grid (var_est =x, r1 = y))
melted = melt(interpolated)
melted$var_est = x
melted$r1 = rep(y, times=rep(length(x), length(y)))
ggplot() + geom_tile(data=melted, aes(x=log10(var_est), y=r1, fill=value)) + scale_fill_gradientn(colours=rev(brewer.pal(9,"RdPu")), limits=c(0, 0.5)) + theme_minimal() + ylim(min(dat2$r1), max(dat2$r1)) + geom_point(data=dat2, aes(x=log10(var_est), y=r1), size= 1.7) + geom_point(data=dat2, aes(x=log10(var_est), y=r1, color=fix), size=1.3) + scale_color_gradientn(colours=rev(brewer.pal(9,"RdPu")), limits=c(0, 0.5))+ theme(panel.grid = element_blank(), panel.border = element_blank()) + scale_x_continuous(breaks = c(unique(log10(dat2$var_est))), labels=formatC(c(unique((dat2$var_est))), format = "e", digits=0), limits=c(-6, -1))+ theme(text = element_text(size=20))+ xlab("fitness variance (log scale)") + ylab("valley depth r1")
```

## Supplementary Fig. 9

```{r include=FALSE}
filepath = "./data_r1_additive/"
first = c(0.9, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1)
second = c(1,5,10,50,100)
dat = extract_data(filepath, first, second, 10000)
```

```{r echo=FALSE}
plots = plot_fix_tunnel(dat$fix_dat, dat$simple_dat, dat$her_dat)
# Fig. 9a
plots$fix

# Fig. 9b
plots$tunnel
```